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Abstract. Using dynamical mean-field theory and the non-crossing approximation as impurity 
solver, we study the response of a Mott insulator to strong dc electric fields. The breakdown 
of the Mott insulating state is triggered by field-induced creation of doublon-hole pairs. In 
a previous investigation, Ref. [1], it was found that the system approaches a long-lived quasi- 
steady state in which the current is time-independent although the number of carriers constantly 
increases. Here we investigate and clarify the nature of this state, which exists only because 
thermalization is slow in the Hubbard model at strong coupling. The current is time-independent 
because doublons and holes have an infinite temperature distribution. Evidence for this fact is 
obtained from spectral functions and by comparing the electric current with the field-induced 
doublon-hole creation rate. Implications to real experiments, in systems with energy dissipation, 
are discussed. 



1. Introduction 

The study of nonequilibrium phenomena in complex strongly correlated systems may reveal 
unexpected physics, which could eventually lead to new ways of tuning and controlling material 
properties on ultrafast timescales. The Mott insulator, a system with partially filled bands in 
which electrons are localized due to the strong Coulomb repulsion, is an ideal example of a 
correlated many-body state, and Mott insulators are common among transition metal oxides 
or organic charge-transfer salts [2]. A possible way to drive these systems out of equilibrium 
is to expose them to a strong (dc) electric field F. The resulting dc response, which is non- 
perturbative in the field and thus cannot be explained by response functions of the equilibrium 
state, is known as the dielectric breakdown. A considerable amount of work in the previous few 
years has been devoted to understand this conceptually simple, and yet theoretically challenging 
phenomenon. Experimentally, nonlinear transport in correlated insulators has been studied 
in both oxides [3, 4] and in organic materials [5, 6]. One observes a strong non-linearity 
in the current-field (j-F) characteristics, with a negative differential resistivity between weak 
current and large current regimes. The "strong-field physics" of correlated materials can also be 
addressed in experiments with ultra-cold atoms in optical lattices, e.g., as an elementary probe 
of the Mott phase [7] . 

In an intuitive picture, the Mott insulating ground state decays through the production 
of doublon-hole pairs in the presence of an electric field: In a scalar potential gauge, the field 



implies a potential energy difference F between two neighboring lattice sites in the field-direction 
(choosing units with lattice spacing a = 1 and electron charge e = 1). In the half-filled Mott 
insulator, with unit occupancy at every site, an electron can thus gain the energy U from the 
field which is needed to form a doublon (doubly occupied site) and leave behind a hole, by 
tunneling over a distance l\j = U/F. In this picture, the tunneling process alone contributes a 
current of the order 



expected to carry a current after they have been created. We will refer to r<jh as the dh-creation 
current. Although its definition may seem ad-hoc here, T^h will turn out to be very useful in 
the analysis of the data below. 

For F <C 1 the tunneling is over many lattice sites and thus occurs at an exponentially small 
rate. An exponential scaling with a threshold field Fth for the electric current, 



and for observables related to the doublon-production rate (with different powers of F in the 
pre-factor) is indeed found quite generically, for the fermionic Hubbard model in one dimension 
[8, 9, 10, 11, 12, 13, 14], in infinite dimensions [1], and also for the bosonic Mott insulator 
[15]. A lot of work aimed at understanding the dependence of F t h on the charge gap A c . 
Motivated by the dielectric breakdown of band insulators [9], which can be explained by the 
Landau-Zener mechanism [16], Oka et al. proposed a description of the Mott breakdown in 
the one-dimensional Hubbard model in terms of Landau-Zener tunneling between many-body 
eigenstates [8, 9, 12, 13, 17]. Both numerical work [8, 9] using exact diagonalization and density- 
matrix renormalization group (DMRG), and an analytical approach involving the ground state 
and the first excited state of the Bethe ansatz solution [12, 13] then suggest a relation Fth ~ A;? 
for small U. Another analytical approach is the solution for the case of one spin- J, electron 

in a spin-f-polarized background, for which one finds F th oc Ac' [14]. DMRG results suggest 
that this behavior holds down to smaller polarization (the Mott insulator being the unpolarized 
case). Furthermore, if the voltage drop is applied over one lattice site rather than linearly, one 
has F t h oc A c [11]. In the infinite-dimensional case [1], F t h clearly increases with A c , but the 
region close to the metal insulator transition (A c — > 0) has so far not been studied. 

In the present paper we do not focus on the value of F t h close to the metal-insulator transition, 
but we want to clarify another fundamental question that arises in connection with studies of the 
dielectric breakdown: Most of the theories mentioned above involve isolated systems. After the 
field is switched on, one observes the emergence of a quasi-steady state in which doublons are 
produced at a more or less time-independent rate, and dc properties of the system are obtained 
from this quasi-steady state. When the electric current is computed independently, one finds 
that j is also time-independent [1, 10, 18]. This seems quite peculiar if doublons and holes, 
whose number is constantly increasing, are interpreted as "charge carriers". The main purpose 
of the present paper is to clarify the nature of this quasi-steady state: We demonstrate that its 
existence is related to a failure of the quasi-equilibrium description, i.e., a lack of thermalization 
(Sec. 3.2). We furthermore clarify the relation of the electric current and the doublon production 
rate, and provide an explanation why the doublons and holes which are generated by the field 
appear to be immobile and do not contribute to the current (Sec. 5). This also leads to a 
better understanding of how temperature (Sec. 5) or generic energy dissipation mechanisms 
(Sec. 6) influence the dielectric breakdown, which is essential for understanding (and designing) 
experiments that can probe the non-perturbative effects of Eq. (2) (Sec. 7). Our analysis builds 
on the results of Ref. [1], in which the dielectric breakdown was studied in the limit of large 
dimensions [19], using dynamical mean-field theory (DMFT) [20]. 





(2) 



2. Model and methods 



Throughout this paper we study the paramagnetic Mott insulating phase in the half-filled 
Hubbard model on a ci-dimensional cubic lattice, 

H= E^(*) c t^ + ^E(^t-5)(^-3)- ( 3 ) 

(ij)a i 

Here cj (ci a ) is the creation (annihilation) operator for an electron with spin a at lattice site Ri, 
U denotes the Coulomb repulsion, and Vij is the matrix element for hopping between nearest 
neighbor sites i and j. We initially prepare the system in thermal equilibrium at temperature 
T = 1//3, and apply a homogeneous electric field F(t) for time t > 0. For convenience, F 
is always pointing along the body diagonal i) = (1 . . . 1)*, which simplifies the DMFT self- 
consistency (see below) [21]. The field is turned on to a value F within a switching time to with 
a given ramp profile r(x), 

F(t) = fjFr(t/t ), (4) 

k — j cos(7rx) + j cos(-7rx) 3 for < X < 1 

r(x) = < 1 , . (5) 

1 tor x > . 

The smooth turn-on of F damps transient currents, but does not influence the long-time behavior 
[1]. To incorporate the field into Eq. (3) we use a gauge with pure vector potential A(t), i.e., 
F(t) = —dtA(t)/c, for which the hopping matrix elements acquire a Peierls phase, 

(t) = V^e ie{R ^ -*0 . (6) 

(For a discussion of various gauges within a tight-binding model, see Ref. [22].) In the limit 
d = oo [19], with rescaled hopping V^ = t*/2yd, the problem can be solved exactly using DMFT 
[20] in its nonequilibrium (Keldysh) formulation [23, 24]. The rescaled nearest neighbor hopping 
t* is used as the unit of energy, such that the density of states is given by p(e) oc exp(— e 2 ). Time 
and field are measured in units of h/t* and t* /ea, respectively (K = 1, e = 1, a = 1). 

The DMFT single-site problem is solved by means of the self-consistent hybridization 
expansion [25]. The entire set of equations for the geometry considered here has been discussed 
in Ref. [26], and details of the numerical implementation of the Keldysh equations are given in 
Ref. [27]. In the present work we use the lowest order of the strong-coupling solver, or non- 
crossing approximation (NCA) [28], which allows us to address certain issues that require long 
simulation times. In order to assess the validity of this approximation, we compare these results 
in Sec. 3.1 to the results of Ref. [1], which were obtained using the second order, or one-crossing 
approximation (OCA) [29]. 

Observables 

From the DMFT solution we directly evaluate [27] expectation values of the double occupancy 
d = (x Tli^i^ii) and the kinetic energy E kin = J2ka e k+A(t)nka), where the band energy 
£fc is the Fourier transform of V^. The current is given by j = Sfco- V k+A(t)^ha) : with the 
band velocity Vf- = d^Ck- One can then verify the general fact that the total internal energy 

Etot = E kin (t) + Ud(t) (7) 

changes according to the equation 



E tot {t)=j(t)F(t). 



(8) 




Figure 1. (a) Time-dependent current j(t) for U = 5 and f3 = 5 and various electric fields F. (b) 
Double occupancy for the same parameters, (c) Current (j)t, averaged over times 40 < t < 50. 
The curve (r<ih)t (open symbols) shows corresponding data for the doublon current Eq. (1), see 
Sec. 5. Inset: Spectral function for the equilibrium state at U = 5 and j3 = 5. 

Because the self-consistent strong-coupling solver is conserving in the sense of Kadanoff and 
Baym, Eq. (8) is satisfied also for the approximate solution. 

In addition to the static observables, we define a spectral function from the Fourier transform 
of the retarded local Green function G rct (t,t') = -i&(t - t')({c(t), c< (t')}}, 

A(u, t) = Im / ds e ius G Ict {t + s,t). (9) 

Jo 

For a general non-equilibrium situation, A(uj, t) is not necessarily positive. However, when the t- 
dependence of A(oj } t) can be neglected, A(u, t) is positive, and it can be related to photoemission 
and inverse photoemission spectra in the usual way [30]. In a quasi-steady state, as we will 
encounter below, these properties remain as long as the £-dependence of A(u, t) is very slow 
compared to the inverse of the scale Auj on which A(u, t) changes as a function of u. Technically, 
a finite cutoff s max limits the frequency resolution. In analogy to Eq. (9) we will also look at 
the density of occupied states (corresponding to the photoemission spectrum), 

A<{uj,t) = -Im dse iuls G < {t + s,t), (10) 

7T Jo 

which is the Fourier transform of the lesser Green function G < (t,t') = i(c^ (t')c(t)}. The 
occupation function 

N(uj,t) = A < (uj,t)/A(uj,t) (11) 
gives the Fermi function f(u>) in equilibrium. 



3. The dielectric breakdown current 

3.1. The current 

In this section we briefly recapitulate the main results of Rcf. [1]. We use the lowest order 
of the strong-coupling expansion (NCA) to compute several quantities related to the dielectric 
breakdown, and compare the outcome to OCA results from Ref. [1]. Figure la shows the 
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Figure 2. (a) Current (j)t/F (averaged over times 40 < t < 50) on a logarithmic scale as a 
function of 1/F for various values of U (/? = 5). In this plot, black solid lines correspond to fits 
with Eq. (2), in the intervals covered by the lines, (b) The value F^, obtained from the linear fits 
(filled symbols). Open symbols show corresponding results from Ref. [1], which were obtained 
within the OCA. (c) Current (j)t/F (averaged over times 20 < t < 30) on a logarithmic scale, 
for large U and F (see discussion in Sec. 4). Vertical lines indicate resonances at integer U/F. 

current j(t) for various electric fields F as a function of time, for U = 5 and /3 = 5. The 
spectral function for these parameters, with a well-developed gap at the Fermi energy, identifies 
the initial state as a Mott insulator (inset of Fig. lc). The temperature 1//3 = 0.2 is still too low 
to allow for a sizable number of thermally excited carriers, such that the linear dc conductivity 
is indistinguishable from zero on the scale of the plot. Nevertheless, the value of j(t) at long 
times is nonzero and it strongly increases with F. (The large signal for t < 3, which is linear 
in F, is the current associated with the polarization of the insulator.) The highly nonlinear 
j-F characteristics of the Mott insulator is revealed by plotting the long-time average (j)t of 
the current against F (Fig. lc). The threshold-like increase of (j)t/F around F = 0.6 is the 
hallmark of the dielectric breakdown of the Mott insulator. 

Following Ref. [1], we fit the j-F curves with Eq. (2) in order to determine the threshold field 
F t h (Fig. 2). Note that Eq. (2) can only be expected to hold asymptotically for small F. For 
large F, the j-F curves can even behave non- monotonously (see, e.g., F > 1 for U = 5), which 
will be explained in Sec. 4 below. On the other hand, the exponentially small current for F <C 1 
cannot easily be resolved numerically. This restricts the fits with Eq. (2) to an intermediate 
range of F, as indicated by the black solid lines in Fig. 2a. The resulting threshold field Fth, 
which is shown in Fig. 2b, decreases when U is decreased from the insulating regime towards 
the metal-insulator transition (which is only a crossover at (3 = 5). 

In Fig. 2b, we have also included OCA data from Ref. [1] for the same parameters (open 
symbols). One can see that the higher-order corrections of the OCA are clearly important to get 
a quantitatively correct value of -Fth- The NCA threshold curve is shifted to smaller interactions 
compared to data from Ref. [1]. This behavior resembles the difference between OCA and 
NCA in equilibrium, e.g., for the location of the metal-insulator transition line. Apart from 
that, however, we find that the physics obtained with NCA throughout the insulating regime is 
qualitatively similar to the OCA solution (compare, e.g., Fig. 2 with Fig. 4c of Ref. [1]). Further 
discussions in this paper will thus only be based on NCA results. This allows us to study in 
detail several relaxation processes with and without dissipation, and to obtain high-resolution 




Figure 3. (a) Time-dependent current j(t) for U = 3 and (3 = 5 and various electric fields 
F. (b) Double occupancy for the same parameters, (c) Double occupancy d(T), and (d), linear 
response conductivity o~dc(T), in thermal equilibrium at U = 3 and temperature T = 1/(3. 



results for the spectral function. The calculation of the latter requires data at long times that 
would be accessible within the OCA only at large numerical expense. 

3.2. Failure of the quasi- equilibrium description 

Although the current in Fig. la becomes almost stationary at long times, the system is not in 
a true steady state. This is already clear from a very general energy consideration, because 
the total energy of any closed system in an external field F must increase at the rate given by 
Eq. (8). Because the energy per particle is bounded from above in a single-band model, a time- 
independent current cannot last forever. In general, the energy increase involves both kinetic 
and interaction energy. The double-occupancy d(t), which is proportional to the interaction 
energy, is plotted in Fig. lb. We find that d(t) indeed increases almost linearly with time, while 
j(t) does not change substantially (compare also Fig. lb of Ref. [1]). This behavior is observed 
(in a quasi-steady manner) for rather long times, although in principle d cannot exceed the 
value 0.5. As mentioned in the introduction, it is this peculiar quasi-steady state from which 
properties of the dielectric breakdown are usually inferred. In this section we investigate to what 
extent a simple quasi-equilibrium description can account for the observed behavior. 

In many situations, it is a valid assumption that a system rapidly evolves to a new equilibrium 
state at elevated temperature after the energy is increased. In accordance with this, the simplest 
possibility to describe an isolated system in an external field is to assume that all its properties 
can be obtained from an effective equilibrium state, with a time-dependent temperature T e s(t) 
which follows from the equation CvT e g(t) = j(t)F(t). In particular, the current at small fields 
would be given by 

iquasi-eq.(t) = F&6.C (^eff) , (12) 

where o"dc is the linear conductivity which becomes nonzero for T e g > 0. This argument was 
found to apply to various correlated systems [31, 32, 33]. In the Hubbard model at weak- 
coupling, e.g., the quasi-equilibrium description works if U exceeds some critical interaction, 
while for smaller interaction scattering of the particles is too slow to establish the equilibrium 
state, and the system performs long-lived Bloch oscillations [33]. 



From the time-independence of j(t) in Fig. la we can already infer that a simple quasi- 
equilibrium description is not valid for the dielectric breakdown at U = 5: Because the linear- 
response current in the insulator depends exponentially on temperature, a quasi-equilibrium 
state which is consistent with the increase of d(t) would imply a current which is far too large. 
In the case of the Mott insulator at U = 5 and F = 1, e.g., the double occupancy increases 
by 0.006 from t = 10 to t = 35 (Fig. lb), while the conductance remains at j/F ~ 0.001. On 
the other hand, in order to increase d by 0.006 in equilibrium, one would have to increase the 
temperature from (3 = 5 (d ~ 0.010) to j3 = 2 [d ph 0.016), but at (3 = 2, the linear response 
conductivity is already more than one order of magnitude larger than the quasi-steady value at 
F = 1 [a dc {/3 = 2) » 0.0155]. 

The lack of thermalization in the Hubbard model at U t* has been encountered previously 
[34, 35], and it implies the existence of various interesting metastable states in the Hubbard 
model [36, 37, 38]. Slow thermalization is ultimately related to the fact that the recombination 
time of doublons and holes is exponentially long for U 3> t* [39, 40], and hence kinetic energy 
and interaction energy cannot efficiently be redistributed. The long lifetime of doublons has been 
suggested earlier as a prerequisite for the existence of a quasi-steady current [1, 18]. Because 
the thermalization time r t h of the double occupancy in the paramagnetic Mott insulator has 
recently been determined within DMFT (cf. Fig. 2b of Ref. [26], for a final temperature (3 = 2), 
we can make the argument more quantitative in the following. For U = 5, a value r t h ~ 4000 
has been found (within NCA), which clearly exceeds the scale of Fig. 1. On the other hand, 
Tth can become as short as a few inverse hoppings when U is in the metal-insulator crossover 
regime, and it will thus be interesting to see what a more rapid thermalization implies for the 
behavior of the current. Figure 3 shows j(t) for U = 3 (where Ref. [26] gives r t h ~ 40 at (3 = 2.) 
In this case the system does indeed no longer establish a time-independent current. Instead, 
j(t) increases with time as long as d(t) is still small (F = 0.1, 0.2, 0.4), and it decreases when 
d(t) approaches the value 1/4 (F = 1.2, 1.4). Qualitatively, this is consistent with the quasi- 
equilibrium behavior: At low temperatures, U = 3 shows an insulating character, such that 
both the double occupancy and the dc-conductivity increase with temperature (Fig. 3c and d). 
At high temperatures, a<± c always decreases with T {a dc — > for T — > oo), while d increases to 
the uncorrelated value d = 1/4. Quantitatively, the current j(t) in the driven system at U = 3 
remains below the quasi-equilibrium value, in particular for large F, where the energy increase 
is faster than the thermalization rate. 

In summary, the above discussion confirms that the existence of a quasi-steady current is 
tied to the lack of thermalization at large U 3> t* . Since the thermalization time depends 
exponentially on U/t*, the loss of the quasi-steady current occurs quite abruptly when U 
is decreased. (A different regime of fast thermalization to infinite temperature has been 
encountered at U = F, when doublon-hole production requires only a single hopping process 
[18].) That there is no steady current at smaller U also explains why we could so far not really 
determine F t ^ down to the metal-insulator transition. From now on, we will focus on regimes of 
large U where the steady current does exist. 

4. Spectral function 

Further insight into the quasi-steady current-carrying state can be gained from the spectral 
function, Eq. (9). In Fig. 4a and b, we plot A{uj,t) for U = 10. The values of F in Fig. 4a 
and b are well below i*th (which is large for U = 10), such that there is essentially no current 
after the initial polarization of the system, and no change of either energy or double occupancy 
within the numerical accuracy. The state is thus truly stationary for all practical purposes, and 
A(u, t) has no measurable i-dependence for t > 10. The spectrum consists of upper and lower 
Hubbard bands, separated by U. For large F, these Hubbard bands split into isolated peaks 
with a spacing Au = F, which are narrowed with respect to the original Hubbard bands. 




Figure 4. (a) Spectrum A(u), t) for t = 10, U = 10, /3 = 5, and various fields, (b) Same data on 
a logarithmic scale. The low-amplitude noise is due to a finite cutoff s max = 40 of the Fourier 
integral Eq. (9). (c) and (d) Spectrum A(u),t) for t = 10, U = 5, j3 = 5. 



The side-bands resemble the Wannier-Stark ladder for an electron in a tight-binding band 
with applied homogeneous field [41, 22]. However, while all single-particle eigenstates in a tight- 
binding model with applied field are rigorously localized, we note that this is not true for the 
many-body situation of a doublon in the Mott insulator. A single doublon or hole which is 
added to the Mott insulator remains mobile even at large F, because its potential energy can 
be transferred to other particles or spin-excitations [42]. Nevertheless, the Hubbard side bands 
still reflect the localization which remains on short times, i.e., the side-peak at uj = U — nF may 
be interpreted as adding an electron at site i into a many-body Wannier-Stark "resonance" that 
is localized at a different site j, with (Rj — Rj)F /\F\ = (Rj — Ri)r] = n. 

At small fields, Wannier Stark side-bands are not clearly resolved, and we observe only a 
broadening of the spectrum (F < 0.8 in Fig. 4b and d). At smaller U, this field-induced 
broadening eventually leads to a filling- in of the gap (Fig. 4c and d), which indicates the 
possibility of the system to create a doublon-hole pair at no cost of energy. The Mott insulating 
ground state becomes unstable in the field. As soon as the Wannier-Stark side-bands split off, 
A(uj = 0) behaves non-monotonously as a function of F. This goes along with a non-monotonous 
j-F curve, and explains the failure of Eq. (2) at large -F. At large U, where the system can 
tolerate strong enough fields such that the side-bands are well resolved, the non-monotonous 
behavior of the current even reveals clear resonances at integer U/F (Fig. 2c), although in 
this regime the steady current is superimposed with more and more irregular and longer lived 
transient oscillations. We note that the case U/F = 1 is particularly interesting because the 




Figure 5. (a) Current for U = 4.5 and high temperature (3 = 3. Inset: The curves of the main 
plot are fit with an exponential j(t) = Aexp(— at) + B, and the rate a is plotted as a function of 
F. The solid line corresponds to the relation a oc F 2 . (b) Current j(t) (solid line) and doublon 
creation current r^hf^) (dashed line) on a logarithmic scale, for U = 4.5 and indicated values of 
P and F. 



degenerate manifold of states is then described by an effective spin model with a quantum phase 
transition [43], but this physics shall not be addressed here. 

5. Influence of temperature 

We now turn to one of the most important discussions of this paper, and ask why doublons and 
holes that are created by field-induced tunneling do not lead to an increase of the current with 
time. Another question to be addressed here is the influence of the initial temperature on the 
dielectric breakdown, which actually turns out to be closely related: For a hot initial state one 
has thermally excited carriers in addition to those which are induced by the field, and yet we 
will find that they do not influence the current at long times. In Fig. 5a we demonstrate this 
behavior for U = 4.5 and a rather high temperature, (3 = 3. For these parameters, the linear 
response current Fade is actually larger than the field-induced current that is expected from 
Eq. (2). However, j(t) is found to be close to Fa^c only at early times, while it decreases almost 
to zero later on. This behavior was already described in Ref. [1], where we also found that the 
final value is close to the zero-temperature field-induced current, although the decay to that 
value is very slow for small F, which hindered a precise determination of the long-time limit for 
most parameters. 

Because we are now using NCA, it is possible to take a closer look at the long-time behavior. 
To get a deeper understanding, it is useful to separately consider the component of the current 
that is related to the increase of the interaction energy, (^£ , j nt )/i ? . This turns out to be just 
the doublon-creation current that was defined ad-hoc in Eq. (1), {f t E int )/F = (1/ F)j- t Ud{t) = 
r<ih(i)- Remarkably, we find that the long time limit of j(t) is very close to T^h, which by itself 
does not change much over the whole time interval (Fig. 5b) . This shows that the current at long 
times is essentially determined by the doublon production rate, while the produced doublons 
themselves seem not to contribute to the current. This fact is also seen in Fig. lc, where we 
plot the long-time average (r<ih)i as a function of F in addition to the time average {j)t- 

To explain this behavior, one may wonder in the first place how charge carriers in correlated 
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Figure 6. Occupation function A < (ui,t) for U = 4.5, f3 = 3 in the upper Hubbard band, for 
various times. F = 0.2 for (a), and F = 0.3 for (b). The curves are compared to the spectral 
function A(u),t) (rescaled by 0.0045), which is almost time-independent for those parameters 
(bold black line). 



systems behave in the presence of large fields. This is well studied, e.g., in the t J model, where 
all terms changing the doublon or hole numbers have been projected out, and in a strong field 
one is left with the motion of quasiparticles. While noninteracting particles in a tight binding 
model are completely localized and perform Bloch oscillations (see, e.g., Ref. [44]), the motion 
of carriers in a many-body system is a competition between Wannier-Stark localization and 
transfer of potential energy to other degrees of freedom, such as spins in the t J model [45, 42] 
or phonons in the Holstein model [47, 46]. Even in a large field, carriers remain mobile (the 
conductivity typically decreases with a power of F), and the decay of the linear response current 
in Fig. 5a cannot be related to Wannier-Stark localization. 

On the other hand, the conductivity will always vanish at infinite temperature, because then 
potential energy can never be passed to other particles, spins, phonons etc. This motivates the 
explanation that for a low doublon production rate, the system can essentially establish a state of 
infinite temperature for the kinetic energy of doublons and holes, while their number is fixed (in 
contrast, infinite temperature for the interaction energy would imply d = 0.25 as in Sec. 3.2). 
More rigorously, to describe the short-time behavior one may project out all terms from the 
Hamiltonian that change the double occupancy (both field induced and interaction induced). 
The resulting model H^h is a generalized tJ model, containing both doublons and holes in an 
external field F. Within this description, doublons and holes might then show quasi-equilibrium 
behavior with a conductivity 0"^, and reach infinite doublon temperature T^h much faster than 
the change of the doublon number happens via tunneling. For T^h — > oo one again has adh = 0, 
and the only remaining current is the ci/i-creation current. 

Although we do not explicitly derive Hdh (we could not solve this model anyhow), there 
are two predictions of this simple argument which one can directly verify: (i) Quite generally, 
from a high-temperature expansion we expect that when T — > oo one has both £7tot ~ 1/T and 
o"dc ~ 1/T. In a quasi-equilibrium description, j(t) should thus decay exponentially with a rate 
oc F~ 2 for F — > 0, due to Eq. (8) [31, 33]. This behavior is indeed seen in the inset of Fig. 
5a. (ii) To directly observe the "heating" of doublons and holes one can look at the occupation 
function, or photoemission spectrum Eq. (10), which is plotted in Fig. 6 for U = 4.5 and various 
times. The fields, F = 0.2 and 0.3, are chosen such that the total weight of A < (u, t) in the upper 
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Figure 7. (a) Comparison of the current at U = 4.5, /3 = 10 with dissipation (solid lines) 
and without dissipation (dashed lines), (b) Comparison of the doublon-creation current (dashed 
lines) and the full electric current j (solid lines) for U = 4.5, /3 = 10, and various fields F 
(A = 0.4). 



Hubbard band, which is roughly related to the total number of doublons, does not change much 
with time. The shape of A < (uj, t), however, changes considerably. In the initial state and shortly 
after the switch-on of the field, the occupation is concentrated at low frequencies, as expected 
for a thermal state with A < (uj) = A{uj)f{uj). As time increases, A < (ui,t) approaches a curve 
that is proportional to the spectral function itself (bold lines in Fig. 6), corresponding to a flat 
distribution, N(uj) = const., characteristic of a state with zero kinetic energy (high-temperature 
state). 

6. Influence of dissipation 

In a condensed matter system, dissipation of energy to lattice-, spin-, and other degrees of 
freedom is unavoidable, and can happen on sub-picosecond timescales. It is thus essential 
to understand which of the findings of the previous section will be robust against such 
processes. The steady states of dissipative systems with dc driving is itself an interesting area 
of research [49, 48]. Within DMFT, dissipation terms can be included phenomenologically, 
either by coupling a local fermionic bath [50, 49, 51] or a bosonic ("phonon") bath [52] at fixed 
temperature. Here we use the latter approach, since it guarantees particle number conservation 
and dissipates only energy. 

Following Ref. [52], the electronic self-energy in this case is supplemented by a 
bath contribution, which is the lowest order diagram for a Holstein-type electron- 
phonon coupling, Sai ss [G] = \ 2 G(t,t')D(t,t'). here, A measures the coupling strength, 
and D(t,t') is the equilibrium propagator for a boson with energy ojq; D(t,t') = 
— iTr[Tc exp(— i J c dtuob'b)b(t)b^(t')]/Z. The temperature 1//3 of the bath is fixed, such that 
the bath has no memory (ujq = t* in the following). It is important to note that we choose 
ojq = t* <C U for the phonon frequency, to prevent opening a new channel for doublon-hole 
recombination via phonon emission. In general, the precise way of including the dissipation 
should not matter too much, and we chose parameters such that the equilibrium physics remains 
almost unchanged by the presence of dissipation. (In contrast, for a fermionic bath in the wide- 
band limit the equilibrium state would be modified due to mid-gap states that are created in 




Figure 8. (a) d{t) during a single-cycle THz-pulse with amplitude Fq, for U = 5 and f3 = 5. 
Inset: The electric field pulse, Eq. (13), with r = 3.9269. (b) Integrated doublon current, 
(r c jh)av/^o = UAcI/Fq (filled symbols), where Ad is the change of the double occupancy over 
the field cycle in panel (a). The same parametrization as in Fig. 2a is used. For comparison, we 
include the dc curve {j)t/F for U = 5, taken from Fig. 2a (open symbols). Both have the same 
slope, related to F t ^. 



the spectrum of the insulator with dissipation [51].) 

Figure 7a compares the current j at U = 4.5, with and without dissipation. In contrast to the 
isolated system (A = 0), the current in a system with dissipation never becomes stationary, but it 
increases more or less linearly with time (A = 0.4). Based on the results of the previous sections, 
the explanation is straightforward: If the system can dissipate energy, doublons and holes will 
never reach T^h = oo, and the distribution N(u, t) never becomes flat. Hence doublons and holes 
contribute to the current at all times. Since their number is increasing due to the field-induced 
doublon-hole production, the current increases with time. The doublon-creation current T^h, on 
the other hand, is not expected to depend strongly on the coupling to phonons with an energy 
ooq <C U. This is in fact true, as apparent from Fig. 7b, where we compare the currents j and 
r^h for dissipative systems with electron-phonon coupling A = 0.4. The qualitative behavior of 
r,ih is the same with and without dissipation, i.e., the value is time-independent and increases 
exponentially with F. 

7. Implication for experiments 

To design an experiment that can potentially measure the threshold behavior, one should account 
for the coupling of the sample to leads [53, 10]. Doublons and holes can then escape through the 
boundary, such that their density in the sample becomes stationary and a true steady state is 
reached. (Otherwise, the quasi-steady state would be quite short-lived: extrapolating the linear 
increase of d in Fig. lb, e.g., gives a time of a few 1000 until the uncorrelated value d = 0.25 
would be reached, corresponding to a few picoseconds for t* = leV.) Nevertheless, in the dc 
setup heating of the system is quite substantial, and the current will be a balance between non- 
perturbative doublon-hole production and dissipative effects. To avoid these complications, one 
could probe the strong field behavior by using short pulses, as already proposed for dielectrics 
(band insulators) [54]. If the pulse frequency Q -C U is below the Mott gap, e.g., in the THz 
range, the doublon production might be described by the dc results to a good approximation. 



If a weak external bias V is applied to the sample in addition to the strong THz pulse, the 
induced charge is collected at the leads. The current, averaged over many THz pulses, would 
then be proportional to the total number Ad of doublons created per pulse. In Fig. 8, we have 
simulated the outcome of such an experiment (as sketched in the inset). Figure 8a shows the 
time-evolution of the double occupancy during a single-cycle pulse 



centered around to = with frequency = 0.4, and a Gaussian envelope (see inset 

in Fig. 8a). Motivated by our results, we define an averaged doublon creation current 
(Fdh)av = U/FqAcI. As seen in Fig. 8b, the threshold behavior of (r<ih)av turns out to be the 
same as in the dc case (Fig. 2). This seems reasonable if the doublon production is determined 
by the largest field during the cycle. 

Using time-resolved photoemission spectroscopy, one might be able to observe the heating 
effect demonstrated in Fig. 6: In a system with strong external bias, but not strong enough 
to lead to a dielectric breakdown, one may suddenly excite carriers created by photo-doping, 
monitor the evolution of the distribution function, and compare it to the result without applied 
field. 

8. Summary and Conclusion 

In conclusion, we have used DMFT with the NCA-based impurity solver to study the behavior 
of a Mott insulator in strong dc electric fields F. The electric current j is given by a contribution 
r^h which can be associated with the doublon-hole production due to field-induced tunneling 
[Eq. (1)], and a contribution related to the conductance of doublons and holes which are either 
thermally excited or induced by the field. For an isolated system, however, these carriers 
accumulate energy from the field and rapidly reach an "infinite temperature" state with zero 
conductivity, as apparent from a flat occupied density of states. This explains the peculiar steady 
state that has been described earlier, in which the current is time-independent (and given by 
r^h) although the number of carriers constantly increases. In an isolated system, the current 
itself is thus a good measure for the field-induced tunneling. In contrast, when dissipation of 
energy is taken into account, carriers cannot reach infinite temperature. The current now does 
depend on the number of doublons, temperature, and the coupling to the environment, whereas 
the doublon-hole creation current T^h still provides an intrinsic measure of the tunneling. For 
fields F < t*, Tdh increases with F with a threshold behavior, Eq. (2). For larger fields, we 
find that the j-F characteristics reveals peaks at integer values of U/F, which are related to 
signatures of Wannier-Stark localization of carriers in the spectrum. 

The threshold field might be determined experimentally by measuring the current induced 
by a field pulse under weak bias, while the heating of induced doublons should be observable 
with time-resolved photoemission. 
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